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Recent years have witnessed tremendous progress in numerical relativity and an ever improv- 
ing performance of ground-based interferometric gravitational wave detectors. In preparation for 
Advanced LIGO and a new era in gravitational wave astronomy, the numerical relativity and grav- 
itational wave data analysis communities are collaborating to ascertain the most useful role for 
numerical relativity waveforms in the detection and characterization of binary black hole coales- 
cences. In this paper, we explore the detectability of equal mass, merging black hole binaries 
with precessing spins and total mass Mt G [80, 350] , using numerical relativity waveforms and 
template-less search algorithms designed for gravitational wave bursts. In particular, we present a 
systematic study using waveforms produced by the MayaKranc code that are added to colored, 
Gaussian noise and analyzed with the Omega burst search algorithm. Detection efficiency is weighed 
against the orientation of one of the black-hole's spin axes. We find a strong correlation between the 
detection efficiency and the radiated energy and angular momentum, and that the inclusion of the 
I = 2, m = ±1, modes, at a minimum, is necessary to account for the full dynamics of precessing 
systems. 
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I. INTRODUCTION 

The detection of gravitational waves is one of the most 
exciting developments expected for the next decade in 
Physics and Astronomy. Ground-based detectors, such 
as LIGO [l| and Virgo Q , have achieved their initial de- 
sign sensitivities and are making progress toward their 
advanced configurations [3UJ]. According to stellar pop- 
ulation synthesis models Q, Advanced LIGO and Ad- 
vanced Virgo will be sensitive to measurable rates of 
compact binary coalescences, and open a new field of 
gravitational-wave astronomy. 

There are many uncertainties in the amount and form 
of gravitational wave signals that will be first detected. 
We refer to gravitational wave bursts as transient sig- 
natures lasting between a few milliseconds and seconds 
in the detector sensitive band, as could be produced by 
core-collapse supernovae or the late stages of the merger 
of two intermediate-mass black holes, as well as serendip- 
itous sources. Burst signatures are typically searched for 
with template-less algorithms that identify coincident ex- 
cess power in multiple detectors [6, 7]. In particular, this 
paper focuses on the detectability of the mergers of com- 
pact objects with total mass Mr £ [80, 350]M Q , whose 
solutions require simulations of the non-linear Einstein 
field equations of General Relativity. The field of numer- 
ical relativity (NR) has produced solutions to the merger 
of two black holes [8l-f29j that are now being used as sig- 
nals to test the performance of analysis routines (30j and 
to calibrate the analytical and phenomcnological wave- 
form families used as template banks in searches [3ll - l33| . 

For systems with total mass greater than ~ 80M Q , 
the final stages of a binary black hole (BBH) coales- 
cence, merger and ringdown, yield most signal-to-noisc 



ratio (SNR) in ground-based detectors. This is the moti- 
vation for the work presented in this paper: a systematic 
study of how the parameters of the coalescence affect 
the detectability of mergers of BBHs in a burst search. 
We test the detectability of mergers in simulated Gaus- 
sian noise at the initial LIGO design sensitivity using the 
Omega algorithm [34[, used by the LIGO- Virgo collab- 



oration to search for gravitational wave bursts, with a 
fixed, single-detector threshold of SNR > 5.5, chosen to 
match the definition of single-detector "trigger" in com- 
pact binary coalescence searches in real LIGO noise [35[ . 
We are interested in the effects of both numerical and 
physical parameters; in this study, we focus on the ef- 
fect of spin and its orientation as well as the number of 
included harmonics and gravitational-wave cycles. 

This paper is organized as follows. In Section [TTl we 
explain how the waveforms are generated and how the 
Omega pipeline burst search algorithm works. In Sec- 
tion IIII1 we explain how we use the NR waveforms to 
generate the simulated signals, and discuss the analysis 
and postprocessing of the data. Results are presented in 
Section [TV] and conclusions in Section IVl 



II. BACKGROUND 
A. NR Waveforms 

The waveforms used in this work were generated 

using the vacuum BBH numerical relativity code 

MayaKranc, also used in previous studies j2ll l36l - [38j . 
MayaKranc adopts the Baumgarte-Shapiro-Shibata- 
Nakamara formulation with moving puncturegauge con- 
ditions [UGll and the KRANC code generator (39|. Bowen- 
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FIG. 1: Spin configuration used in the simulations for this 
work. The masses and spin magnitudes of the black holes are 
equal, and the orbit is initially quasi-circular. 



York extrinsic curvature is used to solve the momentum 
constraint initially poj . and the Hamiltonian constraint 
is solved using the TWOPUNCTURES spectral solver (HI ]. 
This study uses simulations of equal-mass black holes 
initially located on the x-axis, with momentum deter- 
mined b y th e 3rd-order post-Newtonian angular momen- 
tum [42ll43j ]. according to the geometry in Figure [TJ Most 
simulations have an initial separation of d/M = 6.2, al- 
though two sets of simulations with larger initial sepa- 
rations of d/M = 8 and d/M = 10 are also included. 
Our code units are given by M = mi + mi = 2m, 
where m is the horizon mass of the initial black holes. 
In all simulations, the black holes are spinning, with 
S_ = {— a, 0,0}m and S + = {a sin 8, 0, a cos 8}m. We 
vary 8 from 0° to 360° for a given a/m, and we consider 
three spin values: a/m = {0.2,0.4,0.6}. We hold d/M, 
a/m and the total initial mass of the back holes constant 
as we vary 8. 

The grid structure for each run consists of ten levels 
of refinement provided by CARPET [iij . a mesh refine- 
ment package for CACTUS [45|. We use sixth-order spa- 
tial finite differencing; the outer boundaries are located 
at 317M, the finest resolution is M/77 and the wave- 
forms are extracted at 75M. A few waveforms have been 
generated at resolutions of {M/64, M/77, M/90}, with 
convergence consistent with our sixth-order code. The 
short (long) runs exhibit a phase error on the order of 
5 x 10~ 3 (5 x 10 -2 ) radians and an amplitude error of 
ss 2% 5%). Similar accuracy can be expected in all 
runs we performed. 

To characterize the NR runs, in Figure [2] we plot as 
a function of 8 the ratio E Iad /M ADM , that is the to- 
tal radiated energy (E lad ) normalized by the total mass 
of the system given by the Arnowitt-Descr-Misncr mass 
(M ADM ). The radiated energy varies from approximately 
2.8% to almost 4.4 % of the total ADM mass. The min- 
imum occurs when 8 « 180° due to the way the sys- 
tem is configured initially. By keeping the initial sepa- 
ration fixed, the misalignment of the spin with the or- 
bital angular momentum results in 8 — 180° being the 
shortest run with very little inspiral before merger. The 
dashed line was computed using only the dominate mode, 
I = 2, m = ±2, the solid has all modes up to and includ- 
ing I = 6 and the dashed all the m-modes of £ = 2 
(m = 0, m = ±1, m = ±2). Not including all the £ = 2 
modes in computing E rad clearly underestimates the en- 
ergy radiated for some values of 8. Modes other than 
I = 2, however, contribute little to the total radiated 
energy, comparing solid and dashed. 
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2: Radiated energy, normalized by M ADM , plotted ver- 
sus 6 for the series with d/M = 6.2 and a/m = 0.6. The radi- 
ated energy is computed at an extraction radius of r = 75M. 
The solid line is all I < 6, dashed: I = 2, dotted: dominant 
l = % Iml = 2. 
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FIG. 3: The solid curve is the total magnitude of the radiated 
angular momentum normalized by the initial angular momen- 
tum for each run in the series with d/M = 6.2 and a/m = 0.6. 
The dashed curve is the radiated angular momentum in only 
the z-direction with the same normalization. 



In Figure [3l the solid curve is the total magnitude of 
the radiated angular momentum J rad , normalized by the 
initial angular momentum J lnlt , while the dashed curve is 
the radiated angular momentum in the z-direction J^ ad , 
normalized by J™' (the initial angular momentum was 
not fixed for all simulation runs). When the two curves 
align, the simulations are dominated by radiation in the 
z-direction. As we approach 8 = 250°, the other compo- 
nents of the angular momentum contribute significantly. 
The importance of the numeric value of 8 is only in its 
labeling of the binary dynamics, for example that J™ d is 
at a maximum for 8 = 250°. 
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FIG. 4: The dashed line indicates the design sensitivity of 
LIGO's 4km detectors [J; the continuous line is the spectrum 
of simulated noise used in our study. 



B. Burst Searches 

Unmodeled burst searches look for gravitational wave 
transients in LIGO and Virgo data without assumption 
on the waveform morphology [f| @] • This is accomplished 
by searching for statistically significant excess signal en- 
ergy in the strain time series. The advantage of this type 
of analysis is that, because it does not require knowledge 
of the waveform a priori, it has the potential to detect 
signals missed by other searches. The particular burst 
algorithm used for this study is the Omega pipeline [34| . 
This algorithm first whitens the data, then decomposes 
it into a bank of windowed complex exponentials char- 
acterized by a central time r, frequency 0, and quality 
factor Q. The signal is tiled in the T-<p-Q space, and any 
tiles with a normalized energy above some threshold arc 
recorded as possible events. This search is equivalent to a 
matched filtering against whitened data using a template 
bank of sine-Gaussians. 



III. METHODS 

The NR waveforms were exchanged according to the 
data format described in (4(|; the transformation be- 
tween ^4 and strain was done with the method outlined 
in [47| . In our study, NR waveforms are used to simulate 
signals from sources with random sky location and source 
inclination; the same set of simulation is scaled so as to 
effectively produce signals from 27 concentric shells at 
different radial distances from the detector. These simu- 
lated waveform arc then injected into Gaussian noise col- 
ored to mimic the sensitivity of the initial LIGO design, 
shown in Figure 2] The combined signal and noise are 
then fed into Omcg simulated detector output; the 

results are postprocessed to obtain the detection reach 
statistic described below. 

Since this study is meant to provide a relatively broad 




Distance (Mpc) 

FIG. 5: A sample plot of detection efficiency as a function of 
distance to the source. Each point corresponds to the average 
detection efficiency of a spherical shell of radius equal to the 
distance in abscissa. The curve is a sigmoidal fit to the data, 
used to interpolate the distance D 50 % at which the detection 
efficiency falls to 50%, referred to as the reach of the search. 
Note that there is an uncertainty in the value of each point, 
though the error bars are too small to be visible in this plot. 



insight into how efficiently burst searches detect BBH 
coalescences, we allowed some simplifications. First, we 
initially used only the dominant quadrupole modes of 
the waveforms, (£, m) = (2, ±2), and disregarded all oth- 
ers. We found this has a noticeable impact on the wave- 
form detectability; we discuss later the effect of using 
all I = 2 modes. Second, to disentangle waveform effects 
from that of non-stationary noise in real interferometer 
data, we didn't reproduce a full coherent multi-detector 
analysis at fixed false alarm rate, but instead only ana- 
lyzed data from a single, ideal detector, chosen to be one 
of the LIGO 4km detectors, with a fixed SNR detection 
threshold. 

This study is restricted to equal-mass systems with 
a total mass in the 80 — 350 M Q range. Because the 
frequency evolution of the waveforms scales as the to- 
tal mass of the system, this range is chosen so that, in 
the initial detectors, most of the SNR comes from the 
merger and ringdown portions of the waveform. For sys- 
tems with lower total mass, there are enough inspiral 
cycles in the sensitive band to make it appropriate to use 
matched filtering techniques, as opposed to burst search 
algorithms [3(j. For larger values of the total mass, the 
ringdown frequency falls below the 40 Hz high-pass filter 
applied to the data, and the waveform falls below the 
LIGO seismic noise barrier. Our best sensitivity is for 
systems with total mass in the 150 — 200 M Q range. 

Once the Omega algorithms identifies instances of ex- 
cess power in the simulated data, we compare the trigger 
list to the list of injections to identify which were missed 
and found. We use an SNR threshold of 5.5, chosen to 
match the definition of single-detector "trigger" in com- 
pact binary coalescence searches in real LIGO noise [35| . 
This threshold choice is one of convenience, for a study of 
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FIG. 6: Reach for the system discussed in Section IWl where the adimensional spin of each black hole is 0.6, using only the 
(l, m) = (2, ±2) modes of the waveforms. To clarify features at the left and right edges of the plot, four angles were repeated. 
Dotted lines in the contours represent uncertainties from the limited statistics of the runs. Notice that the trend appears to go 
roughly as the total angular momentum of the system: when 8 — 0°, one of the black holes is spin "up", parallel to the orbital 
angular momentum; when 9 = 180° , it is spin "down" , antiparallel to the orbital angular momentum. This trend is broken by 
an asymmetry near 9 — 255° for masses between 125 and 300 Mq; the white rectangle marks the location where this is most 
pronounced. 



the dependence of detectability on various other simula- 
tion parameters. In actual LIGO searches, the detection 
threshold will be determined by the rate of non-Gaussian 
noise transients and what is considered to be an accept- 
able false alarm probability [f| 0]. For each concentric 
shell in our simulated sky, we compute an average de- 
tection efficiency, plot it as a function of the shell radius 
and fit a sigmoid to the data; to quantify the detection 
efficiency of that search, we use this fit to interpolate the 
distance at which the detection efficiency is 50%, herein 
referred to as reach and denoted by D 50 % . In Figure 03 
we show the particular fit for the case a/m = 0.6, 9 = 0°, 
d/M = 6.2. 



IV. RESULTS 

This study focuses on the dependence of detection ef- 
ficiency on the orientation and magnitude of the spins of 
the two merging black holes. 

The system under study consists of two equal-mass 
black holes with the same spin magnitude. In the coordi- 
nate system where the initial orbital angular momentum 
of the system is in the +z direction, the initial spin of 



one black hole is in the — x direction, while the initial 
spin of the second black hole is oriented at an angle 9 
from the +z direction, as sketched in Figure [TJ We ex- 
amine the dependence of detection efficiency on both the 
magnitude of the spins and the parameter 9. 



A. Dependence on total mass 

We begin with a system where the magnitude of the 
adimensional spin of each black hole is 0.6. Numerical 
simulations were generated for 9 ranging from 0° to 345°, 
and the mass was binned into twelve bins, 22.5 Mq wide, 
to produce a contour plot of the reach in the M — 9 plane 
as shown in Figure O For a given value of 9, the reach 
peaks at « 180 Mq. This reflects how the waveform 
frequency structure scales with the system total mass: 
at w 180 Mq, the merger falls into the detector's most 
sensitive frequency band. As the total mass of the system 
either increases or decreases, the waveform portion trig- 
gering Omega correspondingly moves later or earlier into 
the waveform, where less power is emitted, and conse- 
quently the reach is reduced [3(1 EH • This feature is well 
understood, and not unique to this particular system. 
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FIG. 7: Average reach in 80 — 350 M© as a function of 9 
for systems with different initial black hole coordinate sepa- 
rations. The dotted line is d/M = 6.2, dashed is d/M = 8, 
and solid is d/M = 10. 



B. Dependence on spin orientation 

FigurelHlshows that, at a given mass, the reach peaks at 
(9 = 0° and has a minimum near 9 = 180°, though its pre- 
cise location depends on the total mass. Keeping in mind 
that at 9 = 0° one of the black holes is spin "up", par- 
allel to the orbital angular momentum, and at 9 = 180° 
it is spin "down" , antiparallel to the orbital angular mo- 
mentum, the reach appears to go roughly as the radiated 
energy and angular momentum of the system. The reach 
deviates from the total energy radiated in the vicinity of 
9 = 255° and is instead more evocative of J^ ad and of the 
dominant (£,m) = (2, ±2) modes contribution to E rad , 
as seen in Figures [5] and respectively. 



C. Dependence on initial coordinate separation 

Since the waveforms used to generate Figure [5] were 
relatively short, with only a few cycles before merger, we 
repeated the analysis with waveforms of different lengths 
and compared simulations with initial coordinate sepa- 
rations d/M of 6.2, 8, and 10; the system in Figure |B]has 
d/M = 6.2. The results, averaged over the entire mass 
range, are shown in Figure [7] The individual spins of the 
black holes, as well as the orbital angular momentum, 
will in general precess about the total angular momen- 
tum of the system. Therefore, two simulations of different 
lengths with the same initial value of 9 will in general not 
represent the same system. We are interested, however, 
not in confirming point-by-point agreement between the 
simulations, but rather agreement in the overall behav- 
ior of the reach for the different simulations lengths. All 
three simulation sets exhibit the same sinusoidal behav- 
ior observed in the shorter simulations and, to varying 
degree, the asymmetry at 9 = 255°. The effect is more 
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(a) 
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(b) 



FIG. 8: Panel (a) reach as a function of 6 for the system 
sketched in Figure Q] with a/m — 0.2, 0.4, and 0.6, averaged 
over 80 - 350 M : a/m = 0.2 is the dotted line, a/m = 0.4 
the dashed, a/m = 0.6 the solid. The three systems exhibit a 
similar behavior, proportional to the spin magnitude. Panel 
(b) the initial angular momentum is plotted versus 9 for the 



same three spin magnitudes. 



subtle in the simulations with d/M = 8 due to sparse 9 
sampling in the region surrounding the 9 = 255°. 



D. Dependence on spin magnitude 

Next, we compare the behavior of the same system 
for adimensional spin magnitudes 0.2, 0.4, and 0.6. For 
a comparative analysis, we forgo breaking up the data 
into mass bins, and instead show the average reach in 
the 80 — 350 M mass range in Figure 8(a) All three 
systems exhibit the same bulk behavior, with a peak at 
9 = 0° and a minimum at 9 = 180°. The asymmetry in 9 
is evident for spin 0.6, but it does not show in the lower 
spin cases in Figure |8(a) due to the combination of two 
factors: for practical reasons, 9 was sampled less at spin 
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FIG. 9: Average reach in the 80 — 350 Mq mass range, as FIG. 11: Percentage of energy radiated by the non-dominant 
a function of 9, compared to the normalized z-component of modes as a function of 9. 
the radiated angular momentum. Only the dominant m — ±2 
modes are included. 

E. Correlation with radiated energy and angular 
momentum 
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FIG. 10: Reach as a function of 9, restricted to 102.5 — 
125.0 Mq, compared to the normalized ^-component of the 
radiated angular momentum. Only the dominant m = ±2 
modes are included. 



0.4 and 0.2 than at spin 0.6, and the amount of angu- 
lar momentum radiated in the x-dircction decreases with 
decreasing spin, as the precession is decreasing. Overall, 
we find that the variation of the reach is proportional to 
the magnitude of the black hole spins, which is consis- 
tent with the observation that the reach scales with the 
total angular momentum of the system. For comparison, 
Figure |8(b)| shows how the initial angular momentum of 
the system depends on the spin magnitude, as a function 
of the angle 9. 



The asymmetry of the reach versus 9, seen near 9 = 
255° is over-emphasized, since our original analysis did 
not include all the t = 2 modes. Figures 1^1 and [TU1 show 
the reach calculated with only the I = |m| = 2 modes 
versus 9 for two choices of total mass for the system with 
a/m = 0.6 and d/M = 6.2. In Figure [9j the reach is 
averaged over the entire mass range of 80 — 35OM0 while 
in Figure [TU] one bin of total mass is selected: 102.5 — 
125.OM . In each of these plots, JJ ad /4 nit is included 
for reference. 

Comparing the reach calculated from just the £ = 
\m\ = 2 modes with that calculated from all the £ — 2 
modes in Figures and [T31 it is clear that the reach 
asymmetry is alleviated with the inclusion of all the 1 = 1 
modes. 

The dependence of the reach on z-component J^ ad , 
rather than the total radiated angular momentum J rad , 
when only the £ = 2, m = ±2 modes were used in 
the calculations, is explained by Equations (3.22)-(3.24) 
from [4t|, where the radiated momentum in the x and 
y— directions depends on non-zero m =^ £ harmonics. 
Likewise, we get the same effect if one ignores the m ^ I 
modes in calculating the radiated energy. Figure [TT1 gives 
a clear indication of why the value 9 = 255° appears "spe- 
cial" in our choice of initial data by plotting the percent 
energy radiated in the non-dominant £ = 2 modes. 

The reach calculated with all the £ = 2 modes and 
averaged over the 80 — 350M Q mass range, shown in Fig- 
ure 1131 is well represented by the radiated energy. The 
reach, however, does depend on the mass range, as seen 
for instance in Figure [T4l In particular, the lower mass 
bins of data are well described by the radiated angular 
momentum rather than the radiated energy. We specu- 
late that at low masses, the Omega algorithm tends to 
trigger earlier in the waveform, when most of the angular 
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momentum is being radiated. On the other hand, when 
averaged over the full mass-scale, the reach takes on the 
characteristic of radiated energy, shown in Figure [2] as we 
would have expected. 

V. CONCLUSIONS 

We used simulated waveforms for BBH coalescences, 
calculated with the MayaKranc code, to characterize 
the detection efficiency of the LIGO- Virgo Omega burst 
algorithm as a function of the morphological parameters 
of the system and to investigate the algorithm's response 
to bulk features in gravitational waves from BBH merg- 
ers. We chose one particular set of equal-mass, quasi- 
circular, spinning waveforms that allow for precession; 
for fixed initial spin magnitudes and fixed initial separa- 
tion of the black holes in the simulations, we investigated 
the dependence of the detection efficiency on the orien- 
tation of the initial spin vector of one of the black holes, 
parametrized by the angle 9 between such spin and the 
angular momentum of the system. 

As 9 is varied, the initial angular momentum changes 
as in Figure |8(b)| resulting in different dynamics of the 
binary. The detection reach, established with the Omega 
burst algorithm, is sensitive to such differences, peak- 
ing when the initial black-hole spin vector is 9 = 0° and 
dipping near 180°, consistent with the system's initial 
angular momentum. A departure from this trend, noted 
near 9 = 255° , is an indication that non-dominant modes 
are important in these systems. If all the I = 2 modes 
are included in the simulations, the reach is well under- 
stood by the radiated energy when averaging over the 
entire mass range studied, 80 — 350M©; however, when 
the reach is computed at the low-end of the mass range 
its behavior is more evocative of radiated angular mo- 
mentum. We conclude that it is important to account 
for all m ^ I modes to both compute the full radiated 
angular momentum and to measure the detection reach 
of a burst search for precessing BBH systems. 

For non-precessing, aligned-spin, equal-mass systems, 



the higher the magnitude of the black-hole spins, the 
greater the reach, since the system radiates more angu- 
lar momentum [50]. For precessing systems, we find the 
reach depends on the orientation of the spins: higher- 
spin systems have a greater reach when the JJ, mt is at a 
maximum and the spin is aligned with the orbital angular 
momentum, near = 0°, but a lower reach when J™ 1 * is 
at a minimum, when the spin anti- aligned, near 9 = 180°. 
We find the variation of reach with 9 is proportional to 
the magnitude of the black hole spins. For a system with 
two equal magnitude spins of 0.6, the reach variation with 
orientation is ~ 30%, resulting in a ~ 120% difference in 
sensitive volume and in the rate of detectable precessing 
BBH in initial and advanced interferometric detectors. 

The parameter space of BBHs is large and this work 
targets only a portion of it. In future studies, we will 
continue to explore the search algorithm's response to 
the parameters of precessing systems, such as the black 
hole mass ratio, and investigate how degeneracies in the 
detectable parameter space caused by the transition from 
a largely parametrized inspiral phase to a two-parameter 
end state impacts burst algorithms' ability to extract the 
system's physical properties. 
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FIG. 12: Same as Figure [6] but using all £ — 2 modes of the waveforms. Notice that the asymmetry at 9 = 255° is significantly 
reduced by the addition of the additional modes. 
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FIG. 13: Reach as a function of 9 averaged over the entire 
mass range 80—350 Mq for injections using only the dominant 
(£, m) = (2, ±2) modes and injections using all £ = 2 modes. 
Notice that in the vicinity of 9 = 90°, the curves agree very 
well, because this is where the black hole spins are antiparallel 
and cancel, making the contributions of the m — and m = 
±1 modes negligible. As 9 moves away from 90, the additional 
modes become more important, and the curves differ very 
significantly. 



FIG. 14: The solid curve is the reach computed using all £ — 2 
modes, constricted to the lowest mass range 80 — 102.5 Mq. 
This is compared to the radiated angular momentum J lad 
normalized to the initial angular momentum J mlt , the red 
dashed curve. Note that the two curves appear to have very 
similar behavior, indicating that the detection efficiency of 
injection including all £ = 2 modes is roughly proportional to 
the radiated angular momentum. 
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